Introduction
Patterned distributions of individual plant species and species groups are likely to reflect patterns in environmental conditions and constraints imposed on plant survival, growth, and reproductive ability by abiotic and biotic factors, including physiological tolerance of levels of various environmental conditions (e.g., soil moisture, [Ca], wind speed), competition with other plants, and interactions with herbivores and pathogens. If variation in the abundance of focal species characterized by different growth forms or physiologies is correlated with environmental conditions, and not due merely to limited plant dispersal, it can suggest that those forms or physiologies are adapted to different states or levels of those conditions.
Sarracenia is a genus of 11 species of carnivorous plants, occurring on the coastal plains of the eastern United States, that produce pitcher-shaped leaves filled with liquid that function as pitfall traps.
Givnish (personal communication) observed a patchy distribution of S. alata in the extensive wet, sandy pine savannas at Buttercup Flats near Wiggins, Mississippi in the DeSoto National Forest. Such wet pine savannas of the Gulf Coastal plain of the United States contain an extraordinary diversity and abundance of carnivorous plants occurring down slope from mesic longleaf pine (Pinus palustris Mill.) dominated uplands (Brewer 1999). On these sites, S. alata tends to co-occur with crayfish mounds made of clay material (Givnish, pers. obs.), suggesting that the plant’s locally patchy distribution may have originated from a hidden variability in the physical composition of the substrate (Brewer 1999). Specifically, the presence of clay lenses in the subsoil may create localized anoxic zones that prevent nitrogen fixation by soil microbes, leading to wet, N-poor microsites that favor the establishment of carnivorous plants, due to low levels of soil nitrate and hostile, anoxic conditions that can work against root uptake of soil nutrients generally. I therefore aim to test the hypothesis that scattered but hidden clay lenses help generate the patchy distribution of Sarracenia, creating zones that are poor in nitrate and dissolved oxygen.
pattern1
Thirteen transects crossing the described vegetation pattern were set up randomly within in the savanna at Buttercup Flats facing a compass bearing of 160 degrees. Each transect contained a pitcher-absent microsite in the middle and 10-m of pitcher-present microsites on each end. Transects 11, 12, and 13 intersected multiple pitcher-absent patches and were therefore evaluated at more locations compared to the other transects. Soil and tissue nutrients, ground Coverage of all species, and soil clay content clay content were measured at these sites using a stratified random design.
Analysis of Soil Nutrients
Plant Root Simulators (PRS, available from Western Ag) with cation and anion exchange resins be planted in pairs of three in two locations along the transect: 1) a random location at the pitcher-present microsite in the beginning ten meters of the transect and 2) a random location in each pitcher-absent microsite contained within the transect. Probe burial time was 46-47 days. All soil nutrients were measured in μg/cm2.
Each nutrient is evaluated independently for its association with microsite type (pitcher plant-present, pitcher plant-absent).
soilNPK = read.csv(file = "soilNPK.csv", header = TRUE)
soilNPK$Transect = factor(soilNPK$Transect) # interpreted as factor, not numerical values
str(soilNPK)'data.frame': 31 obs. of 6 variables:
$ Microsite: Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
$ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
$ NO3 : num 0.6 5.2 0.9 0.9 3.8 2.8 3.2 1.5 1 1.7 ...
$ NH4 : num 2.3 18.9 1.6 2.6 3.1 1.9 7.3 2 2.8 2.3 ...
$ K : num 128 321 105 130 147 ...
$ P : num 0.9 0.3 0.2 0.3 0.3 0.3 0.4 0.5 0.3 0.9 ...Soil Nitrogren
Soil nitrogen is evaluated as both Total Nitrogen and its components:
Total Nitrogen is treated as the sum of soil nitrate (NO3-) and soil ammonium (NH4+).
Increased clay content would reduce drainage and create local flooding. The resulting anoxic environment would severely limit nitrification and favor denitrification, necesitating alternative nutrient acquisition strategies. We would, therefore, expect lower concentrations of soil nitrogen at clay-rich and (presumably) pitcher-present microsites (Craft and Chiang 2002)
Differences in soil levels of nitrate and ammonium between or within microsites might indicate a preferential aquisition of a given form of nitrogen based on physiological differences between dominant plants (Olde Venterink et al. 2003).
soilNPK$TN <- soilNPK$NO3 + soilNPK$NH4Total Nitrogen
data visualization
layout(matrix(1:2,1,2))
hist(soilNPK$TN, main="", col="tan", breaks=20)
qqnorm(soilNPK$TN, main = '', ylab="soil K")
qqline(soilNPK$TN)
The histogram for total soil nitrogen concentrations is moderately right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles show noticible deviation from theoretical quantiles; Q-Q plot does suggests a significant deviation from normality. Two outliers (large TN values) can be seen in both the histogram and the Q-Q plot. they are on 2 different transects and represent different microsites.
subset(soilNPK, TN>15 ) Microsite Transect NO3 NH4 K P TN
2 - 1 5.2 18.9 321.01 0.3 24.1
29 + 13 0.0 18.5 101.58 1.4 18.5Let’s now just visualize the association between total soil nitrogen and microsite, and also see the pairing by transect (transect 1 is in red, transect 13 is in blue).
transect.color = rep("black",13)
transect.color[c(1,13)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, TN, legend=FALSE, col=transect.color))
plot(TN~Microsite, data=soilNPK)
These two transects show large changes in mean TN (in opposite directions) between microsites (note that Transect 13 has several data points for the “-” microsite).
The outliers are much “less” outliers and the Q-Q plot shows no significant deviation from normality when we consider the log(potassium) values:
layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(TN), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(TN), main = '', ylab="soil log(TN)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(TN), legend=FALSE, col=transect.color))
plot(log(TN)~Microsite, data=soilNPK)
From this, TN values between microsites increase and decrease in an equal number of transects.
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
table(soilNPK[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 2 2
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logTNmeans = with(soilNPK, tapply(log(TN), Microsite:Transect, mean))
TN_minus = logTNmeans[1:13]
TN_plus = logTNmeans[14:26]
TN_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6 -:7
3.1822118 1.2527630 1.5475625 1.2527630 1.3862944 2.0668628 1.0647107
-:8 -:9 -:10 -:11 -:12 -:13
1.3350011 1.6486586 0.8754687 1.7639335 1.2668484 1.6134220 TN_plus +:1 +:2 +:3 +:4 +:5 +:6 +:7
1.0647107 0.9162907 1.9315214 2.3513753 1.3350011 1.5475625 1.7749524
+:8 +:9 +:10 +:11 +:12 +:13
2.5952547 0.7884574 0.2623643 2.1860513 2.6390573 2.9177707 A paired t-test is performed on pitcher plant-absent and pitcher plant-present TN soil concentrations.
t.test(TN_plus, TN_minus, paired = TRUE)
Paired t-test
data: TN_plus and TN_minus
t = 0.55005, df = 12, p-value = 0.5924
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.4678313 0.7838112
sample estimates:
mean of the differences
0.1579899 TN concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.59).
We can do this again with the non-transformed TN values, but there is still no evidence of association (p=0.57).
TNmeans = with(soilNPK, tapply(TN, Microsite:Transect, mean))
TN_minus = TNmeans[1:13]
TN_plus = TNmeans[14:26]
t.test(TN_plus, TN_minus, paired=TRUE)
Paired t-test
data: TN_plus and TN_minus
t = 0.58788, df = 12, p-value = 0.5675
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-3.809498 6.624883
sample estimates:
mean of the differences
1.407692 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with TN (or log(TN)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw TN values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between TN and microsite (p=0.37).
fit <- lm(TN ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: TN
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 244.56 20.380 0.6534 0.7708
Microsite 1 26.81 26.813 0.8596 0.3668
Residuals 17 530.25 31.191 drop1(fit, test="F")Single term deletions
Model:
TN ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 530.25 116.02
Transect 12 254.722 784.98 104.18 0.6805 0.7482
Microsite 1 26.813 557.07 115.55 0.8596 0.3668layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(TN) values this time. Still no evidence of an association (p=0.39). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of TN (untransformed).
fit <- lm(log(TN) ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: log(TN)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 5.0649 0.42207 0.9186 0.5498
Microsite 1 0.3587 0.35873 0.7807 0.3892
Residuals 17 7.8111 0.45948 drop1(fit, test="F")Single term deletions
Model:
log(TN) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 7.8111 -14.732
Transect 12 5.2823 13.0934 -22.718 0.9580 0.5191
Microsite 1 0.3587 8.1699 -15.340 0.7807 0.3892plot(fit)
Ammonium
data visualization
We repeat the same analysis using ammonium (NH4+) the cationic componenet of Total Nitrogen.
layout(matrix(1:2,1,2))
hist(soilNPK$NH4, main="", col="tan", breaks=20)
qqnorm(soilNPK$NH4, main = '', ylab="soil NH4")
qqline(soilNPK$NH4)
The histogran for soil ammonium concentrations is noticibly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not noticible deviation from theoretical quantiles; the Q-Q plot suggests a significant deviation from normality. Four outliers (large ammonium values) can be seen in the Q-Q plots. they are on 4 different transects, representing both microsites equally.
subset(soilNPK, NH4>10 ) Microsite Transect NO3 NH4 K P TN
2 - 1 5.2 18.9 321.01 0.3 24.1
25 - 11 0.0 12.2 110.25 0.8 12.2
26 + 12 0.0 14.0 111.88 0.8 14.0
29 + 13 0.0 18.5 101.58 1.4 18.5Let’s now just visualize the association between ammonium and microsite, and also see the pairing by transect (transect 1 is in red, transects 12 and 13 are in blue).
transect.color = rep("black",13)
transect.color[c(1,13, 12)] = c("red", "blue", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, NH4, legend=FALSE, col=transect.color))
plot(NH4~Microsite, data=soilNPK)
All transects, with the exeptions of the highlighted transects 1, 12, and 13, do not show large differences in ammonium. Ammonium concentrations tend to be slightly greater on the pitcher-present microsite of each transect. Q-Q plots show much less deviation from normality when we consider the log(ammonium) values:
layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(NH4), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(NH4), main = '', ylab="soil log(NH4)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(NH4), legend=FALSE, col=transect.color))
plot(log(NH4)~Microsite, data=soilNPK)
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
table(soilNPK[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 2 2
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logNH4means = with(soilNPK, tapply(log(NH4), Microsite:Transect, mean))
NH4_minus = logNH4means[1:13]
NH4_plus = logNH4means[14:26]
NH4_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6 -:7
2.9391619 0.9555114 0.6418539 0.6931472 0.8329091 1.3609766 0.4054651
-:8 -:9 -:10 -:11 -:12 -:13
0.2623643 0.9555114 0.3364722 1.6098012 0.6404669 1.2644618 NH4_plus +:1 +:2 +:3 +:4 +:5 +:6 +:7
0.8329091 0.4700036 1.1314021 1.9878743 1.0296194 1.3862944 0.6931472
+:8 +:9 +:10 +:11 +:12 +:13
1.5040774 0.5306283 0.2623643 1.8082888 2.6390573 2.9177707 A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.
t.test(NH4_plus, NH4_minus, paired = TRUE)
Paired t-test
data: NH4_plus and NH4_minus
t = 1.1115, df = 12, p-value = 0.2881
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.3172923 0.9781129
sample estimates:
mean of the differences
0.3304103 Ammonium concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.29).
We can do this again with the non-transformed Ammonium values, but there is still no evidence of association (p=0.48).
NH4means = with(soilNPK, tapply(NH4, Microsite:Transect, mean))
NH4_minus = NH4means[1:13]
NH4_plus = NH4means[14:26]
t.test(NH4_plus, NH4_minus, paired=TRUE)
Paired t-test
data: NH4_plus and NH4_minus
t = 0.72334, df = 12, p-value = 0.4833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.991122 5.964199
sample estimates:
mean of the differences
1.486538 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with ammonium (or log(ammonium)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw ammonium values, not transformed. The residuals look fairly normally distributed, but the first residual plot and the Q-Q plot do not look great.
No evidence for association between ammonium and microsite (p=0.31).
fit <- lm(NH4 ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: NH4
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 219.84 18.320 0.7188 0.7159
Microsite 1 28.24 28.241 1.1081 0.3072
Residuals 17 433.27 25.487 drop1(fit, test="F")Single term deletions
Model:
NH4 ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 433.27 109.759
Transect 12 236.380 669.65 99.256 0.7729 0.6700
Microsite 1 28.241 461.52 109.716 1.1081 0.3072layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(ammonium) values this time. Still no evidence of an association (p=0.17). The residual plots look better (more homogeneous variance) and the Q-Q plot show less deviation from normality, so this analysis is preferred over the analysis of ammonium (untransformed).
fit <- lm(log(NH4) ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: log(NH4)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 7.1205 0.59338 1.1050 0.4147
Microsite 1 1.1116 1.11155 2.0699 0.1684
Residuals 17 9.1290 0.53700 drop1(fit, test="F")Single term deletions
Model:
log(NH4) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 9.129 -9.8983
Transect 12 7.8262 16.955 -14.7057 1.2145 0.3480
Microsite 1 1.1116 10.241 -8.3364 2.0699 0.1684plot(fit)
Nitrate
We, again, repeat the same analysis, this time using nitrate (NO3-) the anionic componenet of Total Nitrogen.
data visualization
layout(matrix(1:2,1,2))
hist(soilNPK$NO3, main="", col="tan", breaks=20)
qqnorm(soilNPK$NO3, main = '', ylab="soil NO3")
qqline(soilNPK$NO3)
The histogran for soil nitrate concentrations is moderately right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality (though the first few quantiles look a little weird). Two outliers (large nitrate values) can be seen in the Q-Q plots. they are on 2 different transects and represent different microsites.
subset(soilNPK, NO3>4 ) Microsite Transect NO3 NH4 K P TN
2 - 1 5.2 18.9 321.01 0.3 24.1
15 + 8 8.9 4.5 83.52 0.7 13.4Let’s now just visualize the association between nitrate and microsite, and also see the pairing by transect (transect 1 is in red, transect 8 is in blue).
transect.color = rep("black",13)
transect.color[c(1,8)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, NO3, legend=FALSE, col=transect.color))
plot(NO3~Microsite, data=soilNPK)
Two large changes in nitrate concentrates (those of the highlighted transects) are visible, going in different directions, between microsites.
The outliers are much “less” outliers when we consider the squareroot(nitrate) values. A square root tranformation was chosen for nitrate due to the presence of several zeroes in the data. The Q-Q plot has lost its outliers but still does not look great:
layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(sqrt(NO3), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(sqrt(NO3), main = '', ylab="soil log(NO3)"))
with(soilNPK, interaction.plot(Microsite, Transect, sqrt(NO3), legend=FALSE, col=transect.color))
plot(sqrt(NO3)~Microsite, data=soilNPK)
From this, most transects show slight changes in nitrate concentration across microsites. These changes go in different directions in an equal frequency.
statistical analysis
If we analyze the squart root values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
table(soilNPK[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 2 2
+ 1 1 1 1 1 1 1 1 1 1 1 1 1sqrtNO3means = with(soilNPK, tapply(sqrt(NO3), Microsite:Transect, mean))
NO3_minus = sqrtNO3means[1:13]
NO3_plus = sqrtNO3means[14:26]
NO3_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6 -:7
2.2803509 0.9486833 1.6733201 1.2247449 1.3038405 2.0000000 1.1832160
-:8 -:9 -:10 -:11 -:12 -:13
1.5811388 1.6124515 1.0000000 0.5049690 1.2831928 1.2159615 NO3_plus +:1 +:2 +:3 +:4 +:5 +:6 +:7
0.7745967 0.9486833 1.9493589 1.7888544 1.0000000 0.8366600 1.9748418
+:8 +:9 +:10 +:11 +:12 +:13
2.9832868 0.7071068 0.0000000 1.6733201 0.0000000 0.0000000 A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.
t.test(NO3_plus, NO3_minus, paired = TRUE)
Paired t-test
data: NO3_plus and NO3_minus
t = -0.87215, df = 12, p-value = 0.4002
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.8544107 0.3659245
sample estimates:
mean of the differences
-0.2442431 K concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.92).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.92).
NO3means = with(soilNPK, tapply(NO3, Microsite:Transect, mean))
NO3_minus = NO3means[1:13]
NO3_plus = NO3means[14:26]
t.test(NO3_plus, NO3_minus, paired=TRUE)
Paired t-test
data: NO3_plus and NO3_minus
t = -0.09884, df = 12, p-value = 0.9229
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.816912 1.659220
sample estimates:
mean of the differences
-0.07884615 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with nitrate (or squareroot(nitrte)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw nitrate values, not transformed. The residuals look fairly normally distributed, and the Q-Q plot suggests normality.
No evidence for association between nitrate and microsite (p=0.28).
fit <- lm(NO3 ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: NO3
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 51.902 4.3252 1.3541 0.2765
Microsite 1 0.019 0.0185 0.0058 0.9402
Residuals 17 54.301 3.1942 drop1(fit, test="F")Single term deletions
Model:
NO3 ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 54.301 45.377
Transect 12 51.487 105.788 42.051 1.3432 0.2815
Microsite 1 0.019 54.320 43.388 0.0058 0.9402layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using squareroot(nitrate) values this time. Still no evidence of an association (p=0.34). The residual plots look (slightly) better (more homogeneous variance), so this analysis is preferred over the analysis of nitrate (untransformed).
fit <- lm(sqrt(NO3) ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: sqrt(NO3)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 6.8471 0.57059 1.1713 0.3732
Microsite 1 0.3833 0.38333 0.7869 0.3874
Residuals 17 8.2817 0.48716 drop1(fit, test="F")Single term deletions
Model:
sqrt(NO3) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 8.2817 -12.918
Transect 12 7.1737 15.4555 -17.577 1.2271 0.3409
Microsite 1 0.3833 8.6651 -13.515 0.7869 0.3874plot(fit)
Soil Phosphorus
Soil phosphorus, in contrast to nitrogen, is expected to increase with local inundation, which chemically reduces Fe bound P, releasing Phosphate into the soil (Craft and Chiang 2002).
data visualization
layout(matrix(1:2,1,2))
hist(soilNPK$P, main="", col="tan", breaks=20)
qqnorm(soilNPK$P, main = '', ylab="soil P")
qqline(soilNPK$P)
The histogran for soil phosphorus concentrations is slightly left-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles show noticible deviation from theoretical quantiles; the Q-Q plot suggests significant deviations from normality. One outlier (large P value) can be seen in the Q-Q plots. This data point is on transect 11 and in the “+” microsite: pitcher plants present
subset(soilNPK, P>1.4 ) Microsite Transect NO3 NH4 K P TN
11 + 6 0.7 4 103.23 1.6 4.7Let’s now just visualize the association between phosphorus and microsite, and also see the pairing by transect (transect 9 is in red, transects 6 and 13 is in blue).
transect.color = rep("black",13)
transect.color[c(9,13,6)] = c("red", "blue", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, P, legend=FALSE, col=transect.color))
plot(P~Microsite, data=soilNPK)
The transects with the largest changes in soil (6,9,13) phosphorus do not correspond to the outlier (Transect 11). Soil phosphorus, in general, increases at the pitcher-present site.
The outliers are much “less” outliers when we consider the log(phosphorus) values (but the Q-Q plot still does not look great!):
layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(P), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(P), main = '', ylab="soil log(P)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(P), legend=FALSE, col=transect.color))
plot(log(P)~Microsite, data=soilNPK)
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
table(soilNPK[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 2 2
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logPmeans = with(soilNPK, tapply(log(P), Microsite:Transect, mean))
P_minus = logPmeans[1:13]
P_plus = logPmeans[14:26]
P_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6
-1.20397280 -1.20397280 -1.20397280 -0.69314718 -0.10536052 -0.51082562
-:7 -:8 -:9 -:10 -:11 -:12
0.18232156 -1.20397280 0.18232156 0.09531018 -0.18600488 0.13118213
-:13
-1.06013177 P_plus +:1 +:2 +:3 +:4 +:5 +:6
-0.10536052 -1.60943791 -1.20397280 -0.91629073 -1.20397280 0.47000363
+:7 +:8 +:9 +:10 +:11 +:12
-0.22314355 -0.35667494 -0.69314718 0.09531018 -0.69314718 -0.22314355
+:13
0.33647224 A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.
t.test(P_plus, P_minus, paired = TRUE)
Paired t-test
data: P_plus and P_minus
t = 0.15858, df = 12, p-value = 0.8766
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.4446241 0.5144272
sample estimates:
mean of the differences
0.03490159 P concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.88).
We can do this again with the non-transformed P values, but there is still no evidence of association (p=0.85).
Pmeans = with(soilNPK, tapply(P, Microsite:Transect, mean))
P_minus = Pmeans[1:13]
P_plus = Pmeans[14:26]
t.test(P_plus, P_minus, paired=TRUE)
Paired t-test
data: P_plus and P_minus
t = 0.1976, df = 12, p-value = 0.8467
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.3085107 0.3700491
sample estimates:
mean of the differences
0.03076923 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with P (or log(P)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw K values, not transformed. The residuals look fairly normally distributed.
No evidence for association between K and microsite (p=0.84).
fit <- lm(P ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: P
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 2.1714 0.180951 1.2977 0.3036
Microsite 1 0.0058 0.005796 0.0416 0.8409
Residuals 17 2.3705 0.139443 drop1(fit, test="F")Single term deletions
Model:
P ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 2.3705 -51.697
Transect 12 2.1772 4.5477 -55.500 1.3011 0.3019
Microsite 1 0.0058 2.3763 -53.621 0.0416 0.8409layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(P) values this time. Still no evidence of an association (p=0.86). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of P (untransformed).
fit <- lm(log(P) ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: log(P)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 5.8384 0.48653 1.6440 0.1696
Microsite 1 0.0087 0.00865 0.0292 0.8663
Residuals 17 5.0310 0.29594 drop1(fit, test="F")Single term deletions
Model:
log(P) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 5.0310 -28.369
Transect 12 5.8414 10.8724 -28.481 1.6449 0.1694
Microsite 1 0.0087 5.0396 -30.316 0.0292 0.8663plot(fit)
Soil Potassium
Soil Potassium may be less abundant at microsites with lower clay content due a probable increase in nutrient leeching (Inglett, Reddy, and Corstanje 2005).
data visualization
layout(matrix(1:2,1,2))
hist(soilNPK$K, main="", col="tan", breaks=20)
qqnorm(soilNPK$K, main = '', ylab="soil K")
qqline(soilNPK$K)
The histogran for soil potassium concentrations is slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from theoretical quantiles; Q-Q plot does not suggest any significant deviations from normality. Three outliers (large K values) can be seen in the Q-Q plots. they are on 3 different transects, but all in the “-” microsite: pitcher plants absent (note that one of these transects, #13, has another observation for “-” microsite)
subset(soilNPK, K>200 ) Microsite Transect NO3 NH4 K P TN
2 - 1 5.2 18.9 321.01 0.3 24.1
8 - 4 1.5 2.0 229.97 0.5 3.5
31 - 13 2.2 3.8 225.55 0.6 6.0Let’s now just visualize the association between potassium and microsite, and also see the pairing by transect (transect 1 is in red, transect 4 is in blue).
transect.color = rep("black",13)
transect.color[c(1,4)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(soilNPK, interaction.plot(Microsite, Transect, K, legend=FALSE, col=transect.color))
plot(K~Microsite, data=soilNPK)
The outliers are much “less” outliers when we consider the log(potassium) values:
layout(matrix(1:4,2,2, byrow=T))
with(soilNPK, hist(log(K), main="", col="tan", breaks=20))
with(soilNPK, qqnorm(log(K), main = '', ylab="soil log(K)"))
with(soilNPK, interaction.plot(Microsite, Transect, log(K), legend=FALSE, col=transect.color))
plot(log(K)~Microsite, data=soilNPK)
From this, 4 transects have K higher in the pitcher-present area, and 9 transects have K higher in the pitcher-absent area.
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
table(soilNPK[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 2 2
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logKmeans = with(soilNPK, tapply(log(K), Microsite:Transect, mean))
K_minus = logKmeans[1:13]
K_plus = logKmeans[14:26]
K_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6 -:7 -:8
5.771472 4.868534 4.652435 5.437949 4.959482 3.810876 4.453417 4.704925
-:9 -:10 -:11 -:12 -:13
4.279994 4.211387 4.529404 4.798584 4.838847 K_plus +:1 +:2 +:3 +:4 +:5 +:6 +:7 +:8
4.850858 4.654817 4.988935 4.259718 4.596129 4.636960 4.755915 4.425086
+:9 +:10 +:11 +:12 +:13
4.076859 4.863913 4.371218 4.717427 4.620847 A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.
t.test(K_plus, K_minus, paired = TRUE)
Paired t-test
data: K_plus and K_minus
t = -0.74329, df = 12, p-value = 0.4716
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.4531952 0.2226375
sample estimates:
mean of the differences
-0.1152788 K concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.47).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.23).
Kmeans = with(soilNPK, tapply(K, Microsite:Transect, mean))
K_minus = Kmeans[1:13]
K_plus = Kmeans[14:26]
t.test(K_plus, K_minus, paired=TRUE)
Paired t-test
data: K_plus and K_minus
t = -1.2574, df = 12, p-value = 0.2325
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-72.50872 19.44179
sample estimates:
mean of the differences
-26.53346 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with K (or log(K)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw K values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between K and microsite (p=0.22).
fit <- lm(K ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: K
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 39671 3306.0 1.0959 0.4207
Microsite 1 4940 4939.6 1.6375 0.2179
Residuals 17 51282 3016.6 drop1(fit, test="F")Single term deletions
Model:
K ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 51282 257.75
Transect 12 40938 92220 251.94 1.1309 0.3981
Microsite 1 4940 56222 258.60 1.6375 0.2179layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(K) values this time. Still no evidence of an association (p=0.47). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of K (untransformed).
fit <- lm(log(K) ~ Transect + Microsite, data=soilNPK)
anova(fit)Analysis of Variance Table
Response: log(K)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 2.12012 0.17668 0.9749 0.5064
Microsite 1 0.10048 0.10048 0.5544 0.4667
Residuals 17 3.08095 0.18123 drop1(fit, test="F")Single term deletions
Model:
log(K) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 3.0810 -43.571
Transect 12 2.15191 5.2329 -51.150 0.9895 0.4954
Microsite 1 0.10048 3.1814 -44.576 0.5544 0.4667plot(fit)
Analysis of Sedge tissue nutrients
Aboveground biomass of 6 species of sedges in the Rhynchospora genus (R. megalocarpa, R. compressa, R. gracilenta, R. plumosa, R. rariflora, R. harveyi var harveyi) were collected, dried, and analyzed for % N, P, and K content by the University of Wisconsin Soil and Forage Lab, as a measure of the availability of these nutrients to non-carnivorous plants. All species are C3.
All tissue nutrients were measured in % of Aboveground Dry Mass.
tissuenutrients <- read.table("tissueNPK.csv", header = TRUE, sep=",")
tissuenutrients$Transect = factor(tissuenutrients$Transect)
str(tissuenutrients)'data.frame': 29 obs. of 5 variables:
$ Microsite: Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
$ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
$ TN : num 1.55 1.07 0.63 0.7 1.09 0.81 1.01 0.52 1.08 0.55 ...
$ P : num 0.14 0.0869 0.08 0.06 0.0868 ...
$ K : num 0.67 0.59 0.81 0.61 0.52 0.89 0.6 0.63 0.61 0.9 ...Tissue Nitrogen
Tissue Nitrogen was analyzed only as Total Nitrogen.
Lower soil nitrogen concentrations, combined with inhibition of root uptake on anoxic microsites, could likely result in lower plant tissue N on clay-rich microsites (Olde Venterink et al. 2003).
data visualization
layout(matrix(1:2,1,2))
hist(tissuenutrients$TN, main="", col="tan", breaks=20)
qqnorm(tissuenutrients$TN, main = '', ylab="soil TN")
qqline(tissuenutrients$TN)
The histogram for tissue nitrogen concentrations is moderately right-skewed unimodal. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=29). Sample quantiles do not show noticible deviation from theoretical quantiles; the Q-Q plot does not suggest any significant deviations from normality. One outliers (a large TN value) can be seen in the histogram and Q-Q plot. This data point is on Transect 1 and in the “+” microsite: pitcher plants present
subset(tissuenutrients, TN>1.4 ) Microsite Transect TN P K
1 + 1 1.55 0.14 0.67Let’s now just visualize the association between tissue nitrogen and microsite, and also see the pairing by transect (transect 1 is in red, transect 13 is in blue).
transect.color = rep("black",13)
transect.color[c(1,13)] = c("red", "blue")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, TN, legend=FALSE, col=transect.color))
plot(TN~Microsite, data=tissuenutrients)
The outliers are much “less” outliers when we consider the log(potassium) values:
layout(matrix(1:4,2,2, byrow=T))
with(tissuenutrients, hist(log(TN), main="", col="tan", breaks=20))
with(tissuenutrients, qqnorm(log(TN), main = '', ylab="soil log(TN)"))
with(tissuenutrients, interaction.plot(Microsite, Transect, log(TN), legend=FALSE, col=transect.color))
plot(log(TN)~Microsite, data=tissuenutrients)
From this, TN concentrations tend to change between microsites roughly equally in both directons.
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
table(tissuenutrients[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logTNmeans = with(tissuenutrients, tapply(log(TN), Microsite:Transect, mean))
TN_minus = logTNmeans[1:13]
TN_plus = logTNmeans[14:26]
TN_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6
0.06765865 -0.35667494 -0.21072103 -0.65392647 -0.59783700 0.03922071
-:7 -:8 -:9 -:10 -:11 -:12
0.02955880 -0.06187540 -0.02020271 -0.22314355 0.04306491 -0.31471074
-:13
0.15700375 TN_plus +:1 +:2 +:3 +:4 +:5
0.438254931 -0.462035460 0.086177696 0.009950331 0.076961041
+:6 +:7 +:8 +:9 +:10
0.076961041 0.009950331 -0.030459207 0.086177696 0.122217633
+:11 +:12 +:13
-0.020202707 -0.385662481 -0.301105093 A paired t-test is performed on pitcher plant-absent and pitcher plant-present TN tissue concentrations.
t.test(TN_plus, TN_minus, paired = TRUE)
Paired t-test
data: TN_plus and TN_minus
t = 1.5659, df = 12, p-value = 0.1434
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.05449268 0.33291895
sample estimates:
mean of the differences
0.1392131 TN tissue concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.14).
We can do this again with the non-transformed TN tissue values, but there is still no evidence of association (p=0.12).
TNmeans = with(tissuenutrients, tapply(TN, Microsite:Transect, mean))
TN_minus = TNmeans[1:13]
TN_plus = TNmeans[14:26]
t.test(TN_plus, TN_minus, paired=TRUE)
Paired t-test
data: TN_plus and TN_minus
t = 1.6338, df = 12, p-value = 0.1282
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.04234009 0.29618624
sample estimates:
mean of the differences
0.1269231 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with TN (or log(TN)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw TN values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between TN and microsite (p=0.12).
fit <- lm(TN ~ Transect + Microsite, data=tissuenutrients)
anova(fit)Analysis of Variance Table
Response: TN
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.73255 0.061046 1.7249 0.1583
Microsite 1 0.09506 0.095061 2.6861 0.1220
Residuals 15 0.53086 0.035391 drop1(fit, test="F")Single term deletions
Model:
TN ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.53086 -88.016
Transect 12 0.76379 1.29465 -86.163 1.7985 0.1409
Microsite 1 0.09506 0.62592 -85.239 2.6861 0.1220layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(TN) values this time. Still no evidence of an association (p=0.13). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of TN (untransformed).
fit <- lm(log(TN) ~ Transect + Microsite, data=tissuenutrients)
anova(fit)Analysis of Variance Table
Response: log(TN)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.91612 0.076343 1.6971 0.1655
Microsite 1 0.11542 0.115416 2.5657 0.1301
Residuals 15 0.67477 0.044984 drop1(fit, test="F")Single term deletions
Model:
log(TN) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.67477 -81.060
Transect 12 0.95864 1.63341 -79.422 1.7759 0.1460
Microsite 1 0.11542 0.79018 -78.481 2.5657 0.1301plot(fit)
Tissue Phosphorus
Increased soil P avaliablilty due to frequent inundation would likely result in increased tissue P on clay-rich microsites (Craft and Chiang 2002)
data visualization
layout(matrix(1:2,1,2))
hist(tissuenutrients$P, main="", col="tan", breaks=20)
qqnorm(tissuenutrients$P, main = '', ylab="soil P")
qqline(tissuenutrients$P)
The histogram for tissue phosphorus concentrations is slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality. One outlier (large P value) can be seen in the Q-Q plot and histogram. This data point is on Transect 1 in the “+” microsite: pitcher plants absent
subset(tissuenutrients, P>0.12 ) Microsite Transect TN P K
1 + 1 1.55 0.14 0.67Let’s now just visualize the association between potassium and microsite, and also see the pairing by transect (transect 1 is in red, transects 6 an 13 are in blue).
transect.color = rep("black",13)
transect.color[c(1,13,6)] = c("red", "blue", "blue")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, P, legend=FALSE, col=transect.color))
plot(K~Microsite, data=tissuenutrients)
The outliers are much “less” outliers when we consider the log(phosphorus) values:
layout(matrix(1:4,2,2, byrow=T))
with(tissuenutrients, hist(log(P), main="", col="tan", breaks=20))
with(tissuenutrients, qqnorm(log(P), main = '', ylab="soil log(P)"))
with(tissuenutrients, interaction.plot(Microsite, Transect, log(P), legend=FALSE, col=transect.color))
plot(log(P)~Microsite, data=tissuenutrients)
From this, tissue P concentrations generally decrease in the pitcher-present microsite.
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
table(tissuenutrients[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logPmeans = with(tissuenutrients, tapply(log(P), Microsite:Transect, mean))
P_minus = logPmeans[1:13]
P_plus = logPmeans[14:26]
P_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6 -:7
-2.442560 -2.813411 -2.323911 -2.302585 -2.318024 -2.207275 -2.578154
-:8 -:9 -:10 -:11 -:12 -:13
-2.579459 -2.575116 -2.570648 -2.574808 -2.302585 -2.317587 P_plus +:1 +:2 +:3 +:4 +:5 +:6 +:7
-1.966113 -2.525729 -2.444633 -2.315498 -2.562315 -2.564184 -2.441365
+:8 +:9 +:10 +:11 +:12 +:13
-2.449069 -2.445175 -2.574788 -2.569667 -2.659260 -2.723707 A paired t-test is performed on pitcher plant-absent and pitcher plant-present P tissue concentrations.
t.test(P_plus, P_minus, paired = TRUE)
Paired t-test
data: P_plus and P_minus
t = -0.34997, df = 12, p-value = 0.7324
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.1864111 0.1348144
sample estimates:
mean of the differences
-0.02579836 P concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.73).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.78).
Pmeans = with(tissuenutrients, tapply(P, Microsite:Transect, mean))
P_minus = Pmeans[1:13]
P_plus = Pmeans[14:26]
t.test(P_plus, P_minus, paired=TRUE)
Paired t-test
data: P_plus and P_minus
t = -0.27532, df = 12, p-value = 0.7878
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.01663756 0.01290460
sample estimates:
mean of the differences
-0.001866481 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with P (or log(P)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw P values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between P and microsite (p=0.77).
fit <- lm(P ~ Transect + Microsite, data=tissuenutrients)
anova(fit)Analysis of Variance Table
Response: P
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.0033813 2.8178e-04 1.1784 0.3762
Microsite 1 0.0000212 2.1228e-05 0.0888 0.7698
Residuals 15 0.0035869 2.3912e-04 drop1(fit, test="F")Single term deletions
Model:
P ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.0035869 -232.94
Transect 12 0.0034023 0.0069892 -237.59 1.1857 0.3720
Microsite 1 0.0000212 0.0036081 -234.76 0.0888 0.7698layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(P) values this time. Still no evidence of an association (p=0.71). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of P (untransformed).
fit <- lm(log(P) ~ Transect + Microsite, data=tissuenutrients)
anova(fit)Analysis of Variance Table
Response: log(P)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.39191 0.032659 1.1548 0.3901
Microsite 1 0.00406 0.004060 0.1435 0.7101
Residuals 15 0.42421 0.028281 drop1(fit, test="F")Single term deletions
Model:
log(P) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.42421 -94.520
Transect 12 0.39590 0.82011 -99.403 1.1666 0.3831
Microsite 1 0.00406 0.42827 -96.244 0.1435 0.7101plot(fit)
Tissue Potassium
Tissue Potassium would be expected to increase at clay-rich microsites due to the decreased effect of nutrient leeching presumably resulting in greater K uptake.
data visualization
layout(matrix(1:2,1,2))
hist(tissuenutrients$K, main="", col="tan", breaks=20)
qqnorm(tissuenutrients$K, main = '', ylab="tissue K")
qqline(tissuenutrients$K)
The histogram for tissue potassium concentrations is slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality. Four high outliers (large K values) can be seen in the Q-Q plots. they are on all on different transects and represent the two different microsites equally.
subset(tissuenutrients, K>0.7 ) Microsite Transect TN P K
3 + 2 0.63 0.080000 0.81
6 - 3 0.81 0.097890 0.89
10 - 5 0.55 0.098468 0.90
28 + 13 0.74 0.065631 0.77Let’s now just visualize the association between potassium and microsite, and also see the pairing by transect (transects 2 and 13 are in red, transects 3 and 5 are in blue).
transect.color = rep("black",13)
transect.color[c(2,13,3, 5)] = c("red", "red", "blue", "blue")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, K, legend=FALSE, col=transect.color))
plot(K~Microsite, data=tissuenutrients)
The outliers are much “less” outliers when we consider the log(potassium) values. The Q-Q plot also looks better (though still has a weird high tail):
layout(matrix(1:4,2,2, byrow=T))
with(tissuenutrients, hist(log(K), main="", col="tan", breaks=20))
with(tissuenutrients, qqnorm(log(K), main = '', ylab="tissue log(K)"))
with(tissuenutrients, interaction.plot(Microsite, Transect, log(K), legend=FALSE, col=transect.color))
plot(log(K)~Microsite, data=tissuenutrients)
From this, 2 transects (those in red) have noticibly K higher in the pitcher-present area, and 2 transects (those in blue) have K higher in the pitcher-absent area.
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
A new vector is created specifically for the paired t-test.
table(tissuenutrients[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logKmeans = with(tissuenutrients, tapply(log(K), Microsite:Transect, mean))
K_minus = logKmeans[1:13]
K_plus = logKmeans[14:26]
K_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6
-0.5276327 -0.4942963 -0.1165338 -0.4620355 -0.1053605 -0.6161861
-:7 -:8 -:9 -:10 -:11 -:12
-0.4462871 -0.4780358 -0.6539265 -0.5447272 -0.6598031 -0.7985077
-:13
-0.7339692 K_plus +:1 +:2 +:3 +:4 +:5 +:6
-0.4004776 -0.2107210 -0.6539265 -0.5108256 -0.4942963 -0.6348783
+:7 +:8 +:9 +:10 +:11 +:12
-0.4462871 -0.5276327 -0.6348783 -0.6161861 -0.5621189 -0.8675006
+:13
-0.2613648 A paired t-test is performed on pitcher plant-absent and pitcher plant-present K tissue concentrations.
t.test(K_plus, K_minus, paired = TRUE)
Paired t-test
data: K_plus and K_minus
t = -0.19964, df = 12, p-value = 0.8451
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.1684364 0.1401607
sample estimates:
mean of the differences
-0.01413787 K concentrations are not significantly different between pitcher-present and pitcher-absent microsites (p=0.85).
We can do this again with the non-transformed K values, but there is still no evidence of association (p=0.79).
Kmeans = with(tissuenutrients, tapply(K, Microsite:Transect, mean))
K_minus = Kmeans[1:13]
K_plus = Kmeans[14:26]
t.test(K_plus, K_minus, paired=TRUE)
Paired t-test
data: K_plus and K_minus
t = -0.26949, df = 12, p-value = 0.7921
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.11705568 0.09128645
sample estimates:
mean of the differences
-0.01288462 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with K (or log(K)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw K values, not transformed. The residuals look fairly normally distributed, but the first residual plot does not look great.
No evidence for association between K and microsite (p=0.82).
fit <- lm(K ~ Transect + Microsite, data=tissuenutrients)
anova(fit)Analysis of Variance Table
Response: K
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.20410 0.017008 1.4033 0.2643
Microsite 1 0.00068 0.000680 0.0561 0.8160
Residuals 15 0.18180 0.012120 drop1(fit, test="F")Single term deletions
Model:
K ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.18180 -119.09
Transect 12 0.20460 0.38640 -121.23 1.4068 0.2629
Microsite 1 0.00068 0.18248 -120.98 0.0561 0.8160layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(K) values this time. Still no evidence of an association (p=0.89). The residual plots look better (more homogeneous variance), so this analysis is preferred over the analysis of K (untransformed).
fit <- lm(log(K) ~ Transect + Microsite, data=tissuenutrients)
anova(fit)Analysis of Variance Table
Response: log(K)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.53474 0.044561 1.6578 0.1762
Microsite 1 0.00058 0.000576 0.0214 0.8856
Residuals 15 0.40321 0.026880 drop1(fit, test="F")Single term deletions
Model:
log(K) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.40321 -95.992
Transect 12 0.53394 0.93714 -95.534 1.6553 0.1769
Microsite 1 0.00058 0.40378 -97.951 0.0214 0.8856plot(fit)
NPK Limitation Diagram
N:P, N:K, and P:K ratios from tissue nutrient measurements are compared to critical ratios put forth by (Olde Venterink et al. 2003) to determine whether plants are P or P + N limited (N:P > 14.5, K:P >3.4), K or K + N limited (N:K > 3.1, K:P < 3.4), or N limited (N:P < 14.5, N:K < 2.1). The relationship between N:P:K ratios and type of nutrient limitation is illustrated in the subsequent triaxial graph.
library(ggtern)
plot <- ggtern(data = tissuenutrients, aes(x = TN, y = P, z = K)) +
geom_Tisoprop(value=.3226, color='darkgreen') + # T=Top apex: K/(K+N)=.3226 i.e N/K=1/.3226-1=2.1
geom_Tisoprop(value=.2439, color='lightgreen') + # Top apex: K/(K+N)=.2439 i.e N/K=1/.2439-1=3.1
geom_Risoprop(value=.93548, color='darkblue') + # Right apex P/(P+N)=0.145: N/P=1/.06452-1=14.5
geom_Lisoprop(value=.7727273, color='darkred') + # Left apex K/(K+P)=: P/K=1/0.7727273-1=, K/P=1/0.2941=3.4
geom_point(aes(fill = Microsite),
size = 2,
shape = 21,
color = "black") +
ggtitle("NPK limitation diagram") +
labs(fill = "Microsite") +
theme_rgbw() +
theme(legend.position = c(0,1),
legend.justification = c(1, 1))
plot
Critical ratios are represented by vertical lines.
Plants in both pitcher-present and pitcher-absent microsites are most strongly N-limited (?). CA: I would say that they are most strongly P-limited.
About N limitation:
- All sites have N:P higher than 14.5 (blue isoproportion line). (so none are N limited, correct?).
- None of the N:K values are higher than 3.1 (light green line), but some are higher than 2.1 (darkgreen line). (these would definitely not be N limited, correct?)
About P limitation: all points are super far from the top apex (100% P), close to the horizontal line across that apex, where there is 0% P. More specifically:
- all sites have K:P above 3.4 (red line), and
- all sites have N:P above 4.5 (blue line).
The points are overlapping a lot around the line N:K = 2.1, so it’s hard to see if there are mostly pitcher-present or pitcher-absent sites above/below this line. But let’s tabulate: we get that there are only 4 sites with TN:K above 2.1. Of these 4 sites, 3 are in a pitcher-absent patch (2 along the same transect 11, the third in transect 13), and the 4th is in a pitcher-present area. It goes in the expected direction, but 3 versus 1 is not statistically different from a half/half ratio.
with(tissuenutrients, table(TN/K<2.1, Microsite)) Microsite
- +
FALSE 3 1
TRUE 13 12subset(tissuenutrients, TN/K>2.1) Microsite Transect TN P K
1 + 1 1.55 0.140000 0.67
24 - 11 1.12 0.075815 0.50
25 - 11 1.18 0.075889 0.55
29 - 13 1.17 0.098511 0.48transect.color = rep("black",13)
transect.color[c(2,13,11)] = c("red", "blue","purple")
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, TN/K, legend=F, col=transect.color))
plot(TN/K ~ Microsite, data=tissuenutrients)
From this we see that in 10 of the 13 transect, TN:K is higher in the pitcher-present area than in the pitcher-absent patch. Does that run counter to expectations? (in the plot above and those below, transect 2 is in red, transect 13 is in blue, and transect 11 is in purple.)
If Phosphorous is limiting, we could analyze the ratios P/xxx. For example:
layout(matrix(1:4,2,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/K, legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/TN, legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/(K+TN), legend=F, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, P/(P+K+TN), legend=F, col=transect.color))
We if look at the proportion of P, i.e. P/(P+K+TN), we still don’t find any evidence of association with microsite (pitcher plant present/absent).
fit = lm(P/(P+K+TN) ~ Transect + Microsite, data=tissuenutrients)
drop1(fit, test="F")Single term deletions
Model:
P/(P + K + TN) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.00095217 -271.40
Transect 12 0.00156373 0.00251590 -267.22 2.0528 0.09461 .
Microsite 1 0.00019328 0.00114545 -268.04 3.0448 0.10144
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1layout(matrix(1:4,1,4))
plot(fit)
How about N+P limitation?
layout(matrix(1:2,1,2))
with(tissuenutrients, interaction.plot(Microsite, Transect, (TN+P)/K, legend=T, col=transect.color))
with(tissuenutrients, interaction.plot(Microsite, Transect, (TN+P)/(TN+P+K), legend=F, col=transect.color))
fit = lm((TN+P)/(TN+P+K) ~ Transect + Microsite, data=tissuenutrients)
drop1(fit, test="F")Single term deletions
Model:
(TN + P)/(TN + P + K) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.086167 -140.74
Transect 12 0.087773 0.173940 -144.37 1.2733 0.3246
Microsite 1 0.005863 0.092030 -140.84 1.0206 0.3284layout(matrix(1:4,1,4))
plot(fit)
soilNPKphys = read.csv(file = "soilNPKphys.csv", header = TRUE)
soilNPKphys$Transect = factor(soilNPKphys$Transect) # interpreted as factor, not numerical values
str(soilNPKphys)'data.frame': 31 obs. of 9 variables:
$ Microsite : Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
$ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
$ NO3 : num 0.6 5.2 0.9 0.9 3.8 2.8 3.2 1.5 1 1.7 ...
$ NH4 : num 2.3 18.9 1.6 2.6 3.1 1.9 7.3 2 2.8 2.3 ...
$ K : num 128 321 105 130 147 ...
$ P : num 0.9 0.3 0.2 0.3 0.3 0.3 0.4 0.5 0.3 0.9 ...
$ Clay_Depth: num 17.8 13 12.7 10.2 20 ...
$ NumTowers : int 11 5 4 8 6 0 5 0 5 0 ...
$ Height : num -0.182 -0.682 -1.183 -1.683 -0.682 ...tissueNPKphys = read.csv(file = "tissueNPKphys.csv", header = TRUE)
tissueNPKphys$Transect = factor(tissueNPKphys$Transect) # interpreted as factor, not numerical values
str(tissueNPKphys)'data.frame': 29 obs. of 8 variables:
$ Microsite : Factor w/ 2 levels "-","+": 2 1 2 1 2 1 2 1 2 1 ...
$ Transect : Factor w/ 13 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
$ TN : num 1.55 1.07 0.63 0.7 1.09 0.81 1.01 0.52 1.08 0.55 ...
$ K : num 0.14 0.0869 0.08 0.06 0.0868 ...
$ P : num 0.67 0.59 0.81 0.61 0.52 0.89 0.6 0.63 0.61 0.9 ...
$ Clay_Depth: num 17.8 13 12.7 10.2 20 ...
$ NumTowers : int 11 5 4 8 6 0 5 0 5 0 ...
$ Height : num -0.182 -0.682 -1.183 -1.683 -0.682 ...Analysis of Soil Clay Composition
Clay Depth
Clay depth was measured by placing a soil augar with a 23 cm chamber 69 cm into the soil and recording the maximum height of the clay layer. Measurements were recorded at 1) a random location every 15 m along each transect and 2) following the same stratified random schema and coincident with soil and tissue nutrient measurements. Clay depth measurements were not taken concurrent with the PRS soil nutrient samples to prevent the introduction of local drainage. The number of measurements per transect varies with transect length. All distances are relative.
In general, increased clay content at a given microsite is expected to decrease soil and tissue N levels through inhibition of nitrification and reduction in root functionality as a result of local anoxia resulting from poor drainage at these sites. Conversely, P levels are expected to increase with clay content
layout(matrix(1:2,1,2))
hist(soilNPKphys$Clay_Depth, main="", col="tan", breaks=20)
qqnorm(soilNPKphys$Clay_Depth, main = '', ylab="Clay Depth")
qqline(soilNPKphys$Clay_Depth)
The histogram for clay depth concentrations is slightly left-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality (though it has a wierd high tail). No noticible outliers can be seen in the Q-Q plots.
Let’s now just visualize the association between clay depth and microsite, and also see the pairing by transect.
layout(matrix(1:2,1,2))
with(soilNPKphys, interaction.plot(Microsite, Transect, Clay_Depth, legend=FALSE))
plot(Clay_Depth~Microsite, data=soilNPKphys)
The outliers are much “less” outliers when we consider the log(potassium) values:
layout(matrix(1:4,2,2, byrow=T))
with(soilNPKphys, hist(log(Clay_Depth), main="", col="tan", breaks=20))
with(soilNPKphys, qqnorm(log(Clay_Depth), main = '', ylab="soil log(Clay Depth)"))
with(soilNPKphys, interaction.plot(Microsite, Transect, log(Clay_Depth), legend=FALSE))
plot(log(Clay_Depth)~Microsite, data=soilNPKphys)
There are few large changes in clay depth between microsite but clay depth, in general, decreases on pitcher-present microsites.
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
table(soilNPKphys[1:2]) Transect
Microsite 1 2 3 4 5 6 7 8 9 10 11 12 13
- 1 1 1 1 1 1 1 1 1 1 4 2 2
+ 1 1 1 1 1 1 1 1 1 1 1 1 1logClaymeans = with(soilNPKphys, tapply(log(Clay_Depth), Microsite:Transect, mean))
Clay_minus = logClaymeans[1:13]
Clay_plus = logClaymeans[14:26]
Clay_minus # log values, in fact -:1 -:2 -:3 -:4 -:5 -:6 -:7 -:8
2.566295 2.318458 2.541602 3.170211 2.148559 2.723924 3.209431 3.209431
-:9 -:10 -:11 -:12 -:13
3.209431 2.878074 2.962735 2.958735 2.628579 Clay_plus +:1 +:2 +:3 +:4 +:5 +:6 +:7 +:8
2.878074 2.541602 2.995857 3.101218 2.613923 1.491780 2.349230 2.913166
+:9 +:10 +:11 +:12 +:13
2.379083 1.625311 3.011606 2.541602 2.541602 A paired t-test is performed on pitcher plant-absent and pitcher plant-present K soil concentrations.
t.test(Clay_plus, Clay_minus, paired = TRUE)
Paired t-test
data: Clay_plus and Clay_minus
t = -1.6223, df = 12, p-value = 0.1307
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.63827946 0.09344678
sample estimates:
mean of the differences
-0.2724163 Clay depths are not significantly different between pitcher-present and pitcher-absent microsites (p=0.13).
We can do this again with the non-transformed clay depth values, but there is still no evidence of association (p=0.10).
Claymeans = with(soilNPKphys, tapply(Clay_Depth, Microsite:Transect, mean))
Clay_minus = Claymeans[1:13]
Clay_plus = Claymeans[14:26]
t.test(Clay_plus, Clay_minus, paired=TRUE)
Paired t-test
data: Clay_plus and Clay_minus
t = -1.7628, df = 12, p-value = 0.1034
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-8.3417526 0.8805026
sample estimates:
mean of the differences
-3.730625 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with clay (or log(clay)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw clay depth values, not transformed. The residuals look fairly normally distributed, but the first residual plot and Q-Q plot do not look great.
No evidence for association between K and microsite (p=0.09).
fit <- lm(Clay_Depth ~ Transect + Microsite, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: Clay_Depth
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 460.33 38.361 1.3399 0.28308
Microsite 1 91.65 91.653 3.2012 0.09141 .
Residuals 17 486.72 28.631
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(fit, test="F")Single term deletions
Model:
Clay_Depth ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 486.72 113.36
Transect 12 430.55 917.27 109.01 1.2532 0.32666
Microsite 1 91.65 578.37 116.71 3.2012 0.09141 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1layout(matrix(1:4,1,4, byrow=T))
plot(fit)
We can re-do this analysis, using log(clay) values this time. Still no evidence of an association (p=0.11). The residual and Q-Q plots look better (more homogeneous variance), so this analysis is preferred over the analysis of clay (untransformed).
fit <- lm(log(Clay_Depth) ~ Transect + Microsite, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: log(Clay_Depth)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 2.57569 0.21464 1.3251 0.2901
Microsite 1 0.47465 0.47465 2.9303 0.1051
Residuals 17 2.75368 0.16198 drop1(fit, test="F")Single term deletions
Model:
log(Clay_Depth) ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 2.7537 -47.053
Transect 12 2.38829 5.1420 -51.693 1.2287 0.3400
Microsite 1 0.47465 3.2283 -44.123 2.9303 0.1051plot(fit)
Next, we evaluate the effect of clay depth on soil nutrients using a linear model.
The linear model is constructed with
No evidence for association between soil TN and clay depth (p=0.57), soil P and clay depth (p=0.74), or soil K and clay depth (p=0.96).
We can re-do this analysis, using log(
soilNPKphys$TN <- soilNPKphys$NO3 + soilNPKphys$NH4fit <- lm(TN ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: TN
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 244.56 20.380 0.6341 0.7866
Clay_Depth 1 10.68 10.675 0.3321 0.5720
Residuals 17 546.39 32.141 drop1(fit, test="F")Single term deletions
Model:
TN ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 546.39 116.95
Transect 12 254.603 800.99 104.81 0.6601 0.7652
Clay_Depth 1 10.675 557.07 115.55 0.3321 0.5720plot(fit)
fit <- lm(P ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: P
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 2.17141 0.180951 1.3032 0.3008
Clay_Depth 1 0.01585 0.015852 0.1142 0.7396
Residuals 17 2.36048 0.138852 drop1(fit, test="F")Single term deletions
Model:
P ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 2.3605 -51.829
Transect 12 2.18048 4.5410 -55.546 1.3086 0.2981
Clay_Depth 1 0.01585 2.3763 -53.621 0.1142 0.7396plot(fit)
fit <- lm(K ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: K
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 39671 3306.0 0.9998 0.4878
Clay_Depth 1 7 7.2 0.0022 0.9634
Residuals 17 56215 3306.8 drop1(fit, test="F")Single term deletions
Model:
K ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 56215 260.59
Transect 12 39664 95878 253.14 0.9996 0.4880
Clay_Depth 1 7 56222 258.60 0.0022 0.9634plot(fit)
fit <- lm(log(TN) ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: log(TN)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 5.0649 0.42207 0.8801 0.5806
Clay_Depth 1 0.0171 0.01709 0.0356 0.8525
Residuals 17 8.1528 0.47958 drop1(fit, test="F")Single term deletions
Model:
log(TN) ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 8.1528 -13.405
Transect 12 4.7758 12.9286 -23.111 0.8299 0.6220
Clay_Depth 1 0.0171 8.1699 -15.340 0.0356 0.8525plot(fit)
fit <- lm(log(P) ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: log(P)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 5.8384 0.48653 1.6629 0.1643
Clay_Depth 1 0.0656 0.06564 0.2243 0.6418
Residuals 17 4.9740 0.29259 drop1(fit, test="F")Single term deletions
Model:
log(P) ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 4.9740 -28.723
Transect 12 5.8773 10.8513 -28.541 1.6739 0.1612
Clay_Depth 1 0.0656 5.0396 -30.316 0.2243 0.6418plot(fit)
fit <- lm(log(K) ~ Transect + Clay_Depth, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: log(K)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 2.1201 0.176677 0.9447 0.5294
Clay_Depth 1 0.0020 0.001986 0.0106 0.9191
Residuals 17 3.1795 0.187027 drop1(fit, test="F")Single term deletions
Model:
log(K) ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 3.1795 -42.596
Transect 12 2.12152 5.3010 -50.749 0.9453 0.5289
Clay_Depth 1 0.00199 3.1814 -44.576 0.0106 0.9191plot(fit)
We next explore the interaction between tissue nutrient concentrations and clay content using the same linear model design:
The linear model is constructed with clay (or log(clay)) as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw clay depth values, not transformed. No evidence for association between tissue TN and clay depth (p=0.57), tissue P and clay depth (p=0.40), or tissue K and clay depth (p=0.36).
We can re-do this analysis, using log(
fit <- lm(TN ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: TN
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.73255 0.061046 1.4956 0.2282
Clay_Depth 1 0.01365 0.013651 0.3344 0.5716
Residuals 15 0.61227 0.040818 drop1(fit, test="F")Single term deletions
Model:
TN ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.61227 -83.878
Transect 12 0.73140 1.34367 -85.085 1.4932 0.2291
Clay_Depth 1 0.01365 0.62592 -85.239 0.3344 0.5716plot(fit)
fit <- lm(P ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: P
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.204099 0.0170083 1.4681 0.2384
Clay_Depth 1 0.008705 0.0087051 0.7514 0.3997
Residuals 15 0.173775 0.0115850 drop1(fit, test="F")Single term deletions
Model:
P ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.17377 -120.40
Transect 12 0.188142 0.36192 -123.12 1.3533 0.2861
Clay_Depth 1 0.008705 0.18248 -120.98 0.7514 0.3997plot(fit)
fit <- lm(K ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: K
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.0033813 0.00028178 1.2398 0.3421
Clay_Depth 1 0.0001989 0.00019888 0.8750 0.3644
Residuals 15 0.0034092 0.00022728 drop1(fit, test="F")Single term deletions
Model:
K ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.0034092 -234.41
Transect 12 0.0034610 0.0068702 -238.09 1.269 0.3268
Clay_Depth 1 0.0001989 0.0036081 -234.76 0.875 0.3644plot(fit)
fit <- lm(log(TN) ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: log(TN)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.91612 0.076343 1.4766 0.2352
Clay_Depth 1 0.01463 0.014632 0.2830 0.6025
Residuals 15 0.77555 0.051703 drop1(fit, test="F")Single term deletions
Model:
log(TN) ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.77555 -77.023
Transect 12 0.91064 1.68619 -78.500 1.4677 0.2386
Clay_Depth 1 0.01463 0.79018 -78.481 0.2830 0.6025plot(fit)
fit <- lm(log(P) ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: log(P)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.53474 0.044561 1.7235 0.1587
Clay_Depth 1 0.01597 0.015966 0.6175 0.4442
Residuals 15 0.38782 0.025855 drop1(fit, test="F")Single term deletions
Model:
log(P) ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.38782 -97.121
Transect 12 0.50210 0.88992 -97.034 1.6184 0.1876
Clay_Depth 1 0.01597 0.40378 -97.951 0.6175 0.4442plot(fit)
fit <- lm(log(K) ~ Transect + Clay_Depth, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: log(K)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.39191 0.032659 1.2061 0.3605
Clay_Depth 1 0.02208 0.022078 0.8153 0.3808
Residuals 15 0.40619 0.027079 drop1(fit, test="F")Single term deletions
Model:
log(K) ~ Transect + Clay_Depth
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.40619 -95.779
Transect 12 0.39659 0.80278 -100.022 1.2205 0.3525
Clay_Depth 1 0.02208 0.42827 -96.244 0.8153 0.3808plot(fit)
Crayfish Burrow Exits as a proxy for clay content
Crayfish, Fallicambarus spp., abundant in wet pine savannas create distinctive soil changes by constructing “chimney” and “mound” structures as antechambers extensive underground burrows. Both structures (refered to collectively as towers), preferentially constructed based on time of year and weather, consist of silt or clay subsoil and therefore serve as a reasonable approximation of clay content (Brewer 1999).
Crayfish tower abundance was recorded coincident with clay depth by counting the number of towers visible within a 1 X 1 m PVC square directly to the right of the transect.
We explore tower abundance as a proxy and a possibly more accurate measure for soil clay content compared to relative clay depth by applying on a Poisson regression using number of towers as the dependent variable and Transect and Microsite as the first and second factors respectively.
library(MASS)
fit <- glm(soilNPKphys$NumTowers ~ soilNPKphys$Transect + soilNPKphys$Microsite, family=poisson)
summary(fit)
Call:
glm(formula = soilNPKphys$NumTowers ~ soilNPKphys$Transect +
soilNPKphys$Microsite, family = poisson)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.29640 -1.60832 -0.00022 0.86316 1.79748
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.66268 0.27875 5.965 2.45e-09 ***
soilNPKphys$Transect2 -0.28768 0.38188 -0.753 0.451253
soilNPKphys$Transect3 -0.98083 0.47871 -2.049 0.040473 *
soilNPKphys$Transect4 -1.16315 0.51235 -2.270 0.023193 *
soilNPKphys$Transect5 -1.16315 0.51235 -2.270 0.023193 *
soilNPKphys$Transect6 -0.98083 0.47871 -2.049 0.040473 *
soilNPKphys$Transect7 -0.69315 0.43301 -1.601 0.109431
soilNPKphys$Transect8 -0.98083 0.47871 -2.049 0.040473 *
soilNPKphys$Transect9 -0.37469 0.39167 -0.957 0.338747
soilNPKphys$Transect10 0.22314 0.33541 0.665 0.505868
soilNPKphys$Transect11 0.09524 0.30794 0.309 0.757105
soilNPKphys$Transect12 -19.23701 1977.86986 -0.010 0.992240
soilNPKphys$Transect13 -1.95885 0.62989 -3.110 0.001872 **
soilNPKphys$Microsite+ 0.71004 0.18392 3.861 0.000113 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 132.713 on 30 degrees of freedom
Residual deviance: 49.457 on 17 degrees of freedom
AIC: 152.51
Number of Fisher Scoring iterations: 15Number of towers shows significant interaction with microsite.
Height as a possible controller of nutrient limitation
It is possible that since elevated microsites would be subject to greater leeching and would therefore be more nutrient poor compared to microsites lower in the landscape irrespective of soil composition. Height relative to the beginning of the transect was approximated with a range finder according to standard practices. Measurements were taken following the same stratified random schema and coincident as soil and tissue nutrient measurements.
layout(matrix(1:2,1,2))
hist(soilNPKphys$Height, main="", col="tan", breaks=20)
qqnorm(soilNPKphys$Height, main = '', ylab="Height")
qqline(soilNPKphys$Height)
The histogram for height concentrations is very slightly right-skewed. Q-Q plots were used to test for normality because the Shapiro-Wilk test is not powerful at small sample size (n=31). Sample quantiles do not show noticible deviation from Q-Q plot does not suggest any significant deviations from normality. No noticible outliers can be seen in the Q-Q plots.
Let’s now just visualize the association between height and microsite, and also see the pairing by transect.
layout(matrix(1:2,1,2))
with(soilNPKphys, interaction.plot(Microsite, Transect, Height, legend=FALSE))
plot(Height~Microsite, data=soilNPKphys)
statistical analysis
If we analyze the log values, it seems that (at first sight) we can use a t-test or a 2-way ANOVA-based method. We will double-check the normality of residuals after we get these residuals.
At transects 11, 12, and 13, multiple data points (for soil nutrients only) were collected for pitcher plant-absent microsites. It is necessary to average these pitcher plant-absent data points for their respective transects to perform a t-test (which requires a single value per transect per microsite). This averaging is not necessary with a 2-way ANOVA-based approach.
A new vector is created specifically for the paired t-test.
Height is not significantly different between pitcher-present and pitcher-absent microsites (p=0.11).
Heightmeans = with(soilNPKphys, tapply(Height, Microsite:Transect, mean))
Height_minus = Heightmeans[1:13]
Height_plus = Heightmeans[14:26]
t.test(Height_plus, Height_minus, paired=TRUE)
Paired t-test
data: Height_plus and Height_minus
t = 1.7339, df = 12, p-value = 0.1085
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.07649503 0.67264888
sample estimates:
mean of the differences
0.2980769 A linear model is also constructed as both a supplement and a comparison to the paired t-test. We get practically the same p-values, but we also get to see diagnostic plots to confirm that assumptions are met reasonably well.
The linear model is constructed with height as the dependent variable, transect as the first factor, and microsite as the second factor. We first use the raw height values, not transformed.
No evidence for association between K and microsite (p=0.24).
fit <- lm(Height ~ Transect + Microsite, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: Height
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 36.929 3.07740 7.6375 0.0001082 ***
Microsite 1 0.608 0.60845 1.5100 0.2358714
Residuals 17 6.850 0.40293
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(fit, test="F")Single term deletions
Model:
Height ~ Transect + Microsite
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 6.850 -18.802
Transect 12 37.322 44.172 14.977 7.7188 0.000101 ***
Microsite 1 0.608 7.458 -18.164 1.5100 0.235871
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1layout(matrix(1:4,1,4, byrow=T))
plot(fit)
Next, we evaluate the effect of height on soil nutrients using a linear model.
The linear model is constructed with
We first use the raw nutrient values, not transformed.
No evidence for association between soil TN and height (p=0.57), soil P and height (p=0.74), or soil K and height (p=0.96).
We can re-do this analysis, using log(
fit <- lm(TN ~ Transect + Height, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: TN
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 244.56 20.380 0.6489 0.7744
Height 1 23.19 23.185 0.7383 0.4022
Residuals 17 533.88 31.405 drop1(fit, test="F")Single term deletions
Model:
TN ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 533.88 116.23
Transect 12 264.680 798.56 104.71 0.7023 0.7298
Height 1 23.185 557.07 115.55 0.7383 0.4022plot(fit)
fit <- lm(P ~ Transect + Height, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: P
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 2.17141 0.180951 1.3068 0.2990
Height 1 0.02236 0.022356 0.1614 0.6928
Residuals 17 2.35398 0.138469 drop1(fit, test="F")Single term deletions
Model:
P ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 2.3540 -51.914
Transect 12 1.72934 4.0833 -58.839 1.0407 0.4583
Height 1 0.02236 2.3763 -53.621 0.1614 0.6928plot(fit)
fit <- lm(K ~ Transect + Height, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: K
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 39671 3306.0 1.2092 0.35100
Height 1 9742 9742.2 3.5632 0.07627 .
Residuals 17 46480 2734.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(fit, test="F")Single term deletions
Model:
K ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 46480 254.70
Transect 12 43201 89681 251.07 1.3167 0.29415
Height 1 9742 56222 258.60 3.5632 0.07627 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(fit)
fit <- lm(log(TN) ~ Transect + Height, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: log(TN)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 5.0649 0.42207 0.9013 0.5635
Height 1 0.2093 0.20930 0.4470 0.5128
Residuals 17 7.9606 0.46827 drop1(fit, test="F")Single term deletions
Model:
log(TN) ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 7.9606 -14.144
Transect 12 5.2281 13.1886 -22.494 0.9304 0.5405
Height 1 0.2093 8.1699 -15.340 0.4470 0.5128plot(fit)
fit <- lm(log(P) ~ Transect + Height, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: log(P)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 5.8384 0.48653 1.6458 0.1691
Height 1 0.0141 0.01410 0.0477 0.8297
Residuals 17 5.0256 0.29562 drop1(fit, test="F")Single term deletions
Model:
log(P) ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 5.0256 -28.403
Transect 12 4.8323 9.8578 -31.517 1.3622 0.2727
Height 1 0.0141 5.0396 -30.316 0.0477 0.8297plot(fit)
fit <- lm(log(K) ~ Transect + Height, data=soilNPKphys)
anova(fit)Analysis of Variance Table
Response: log(K)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 2.12012 0.17668 1.0747 0.4349
Height 1 0.38671 0.38671 2.3523 0.1435
Residuals 17 2.79473 0.16440 drop1(fit, test="F")Single term deletions
Model:
log(K) ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 2.7947 -46.594
Transect 12 2.12893 4.9237 -53.038 1.0792 0.4319
Height 1 0.38671 3.1814 -44.576 2.3523 0.1435plot(fit)
We next explore the interaction between tissue nutrient concentrations and height using the same linear model design:
The linear model is constructed with
fit <- lm(TN ~ Transect + Height, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: TN
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.73255 0.061046 1.4680 0.2385
Height 1 0.00216 0.002155 0.0518 0.8230
Residuals 15 0.62376 0.041584 drop1(fit, test="F")Single term deletions
Model:
TN ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.62376 -83.339
Transect 12 0.73455 1.35831 -84.771 1.4720 0.2369
Height 1 0.00216 0.62592 -85.239 0.0518 0.8230plot(fit)
fit <- lm(P ~ Transect + Height, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: P
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.20410 0.0170083 1.4780 0.2347
Height 1 0.00987 0.0098698 0.8577 0.3690
Residuals 15 0.17261 0.0115073 drop1(fit, test="F")Single term deletions
Model:
P ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.17261 -120.60
Transect 12 0.20314 0.37575 -122.04 1.4711 0.2373
Height 1 0.00987 0.18248 -120.98 0.8577 0.3690plot(fit)
fit <- lm(K ~ Transect + Height, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: K
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.0033813 2.8178e-04 1.1902 0.3694
Height 1 0.0000570 5.6983e-05 0.2407 0.6308
Residuals 15 0.0035511 2.3674e-04 drop1(fit, test="F")Single term deletions
Model:
K ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.0035511 -233.23
Transect 12 0.0034226 0.0069737 -237.65 1.2048 0.3612
Height 1 0.0000570 0.0036081 -234.76 0.2407 0.6308plot(fit)
fit <- lm(log(TN) ~ Transect + Height, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: log(TN)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.91612 0.076343 1.4520 0.2446
Height 1 0.00149 0.001488 0.0283 0.8686
Residuals 15 0.78869 0.052580 drop1(fit, test="F")Single term deletions
Model:
log(TN) ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.78869 -76.536
Transect 12 0.91695 1.70564 -78.167 1.4533 0.2441
Height 1 0.00149 0.79018 -78.481 0.0283 0.8686plot(fit)
fit <- lm(log(P) ~ Transect + Height, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: log(P)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.53474 0.044561 1.7900 0.1428
Height 1 0.03036 0.030360 1.2195 0.2869
Residuals 15 0.37342 0.024895 drop1(fit, test="F")Single term deletions
Model:
log(P) ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.37342 -98.218
Transect 12 0.52666 0.90008 -96.704 1.7630 0.1490
Height 1 0.03036 0.40378 -97.951 1.2195 0.2869plot(fit)
fit <- lm(log(K) ~ Transect + Height, data=tissueNPKphys)
anova(fit)Analysis of Variance Table
Response: log(K)
Df Sum Sq Mean Sq F value Pr(>F)
Transect 12 0.39191 0.032659 1.1754 0.3780
Height 1 0.01149 0.011494 0.4137 0.5298
Residuals 15 0.41677 0.027785 drop1(fit, test="F")Single term deletions
Model:
log(K) ~ Transect + Height
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.41677 -95.033
Transect 12 0.40141 0.81818 -99.471 1.2039 0.3617
Height 1 0.01149 0.42827 -96.244 0.4137 0.5298plot(fit)
Species Distribution
Species richness and abundance were measured by estimating percent ground cover of each species along each transect 1) at a random location every 15 m along each transect coincident with clay depth measurements and 2) following the same stratified random schema and coincident soil and tissue nutrient measurements.
Species distribution and productivity is influenced by nutrient limitation, with some plants having special adaptations (e.g. symbiotic N2 fixation, root secretion of phosphatases, etc.) that favor nutrient-poor microsites (Olde Venterink et al. 2003).
Clay content, as a proxy for local anoxia, may also influence species distribution according to species flood tolerance.
We use Nonmetric Multidimensional scaling as the ordination method.
First, we look at the species distributions on sites coincident with clay depth data:
speciesClay = read.csv(file = "speciesAbundanceClay.csv", header = TRUE)
library(vegan)Loading required package: permuteLoading required package: latticeThis is vegan 2.4-3ord <- metaMDS(speciesClay[9:35], distance = "jaccard",k=2,trymax=1000,autotransform=TRUE,expand=FALSE, plot=FALSE)Run 0 stress 0.1720377
Run 1 stress 0.1923467
Run 2 stress 0.1730232
Run 3 stress 0.172013
... New best solution
... Procrustes: rmse 0.01030351 max resid 0.07533652
Run 4 stress 0.1988807
Run 5 stress 0.1717517
... New best solution
... Procrustes: rmse 0.02555623 max resid 0.1629871
Run 6 stress 0.1730221
Run 7 stress 0.1921792
Run 8 stress 0.1996325
Run 9 stress 0.1720869
... Procrustes: rmse 0.02679115 max resid 0.1664139
Run 10 stress 0.1730208
Run 11 stress 0.1726898
Run 12 stress 0.1730249
Run 13 stress 0.1851641
Run 14 stress 0.1727169
Run 15 stress 0.1730269
Run 16 stress 0.1720881
... Procrustes: rmse 0.02686078 max resid 0.166446
Run 17 stress 0.1721017
... Procrustes: rmse 0.02744538 max resid 0.1667545
Run 18 stress 0.1729556
Run 19 stress 0.1724165
Run 20 stress 0.1908043
Run 21 stress 0.196815
Run 22 stress 0.1717445
... New best solution
... Procrustes: rmse 0.001226186 max resid 0.006088162
... Similar to previous best
*** Solution reachedscl <- 3 ## scaling == 3
colvec <- c("red2", "mediumblue")
plot(ord, type = "n")
with(speciesClay, points(ord, display = "sites", col = colvec,
pch = 21, bg = colvec))
text(ord, display = "sites", cex = 0.8, col = "darkcyan")
with(speciesClay, legend("topleft", levels(Microsite), bty = "n",
col = colvec, pch = 21, pt.bg = colvec))
adonis(speciesClay[9:35]~speciesClay$Microsite,method="jaccard")
Call:
adonis(formula = speciesClay[9:35] ~ speciesClay$Microsite, method = "jaccard")
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
speciesClay$Microsite 1 2.225 2.22499 7.3231 0.11042 0.001 ***
Residuals 59 17.926 0.30383 0.88958
Total 60 20.151 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(speciesClay[11:31]~speciesClay$Clay_cm,method="jaccard")
Call:
adonis(formula = speciesClay[11:31] ~ speciesClay$Clay_cm, method = "jaccard")
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
speciesClay$Clay_cm 1 0.2966 0.29656 0.89232 0.0149 0.544
Residuals 59 19.6085 0.33235 0.9851
Total 60 19.9051 1.0000 Species show significant grouping by microsite (p=0.001) but not by clay depth (p=0.793)
Next, we look at the species distributions on sites coincident with clay depth data:
speciesNPK = read.csv(file = "speciesAbundanceNPK.csv", header = TRUE)
library(vegan)
ord <- metaMDS(speciesNPK[11:31], distance = "jaccard",k=2,trymax=1000,autotransform=TRUE,expand=FALSE, plot=FALSE)Run 0 stress 0.1658997
Run 1 stress 0.1742643
Run 2 stress 0.1615547
... New best solution
... Procrustes: rmse 0.0727493 max resid 0.2859159
Run 3 stress 0.180824
Run 4 stress 0.1895995
Run 5 stress 0.1606823
... New best solution
... Procrustes: rmse 0.04366153 max resid 0.1825723
Run 6 stress 0.173678
Run 7 stress 0.1673051
Run 8 stress 0.1791267
Run 9 stress 0.1791192
Run 10 stress 0.1592366
... New best solution
... Procrustes: rmse 0.05993994 max resid 0.2076739
Run 11 stress 0.1592368
... Procrustes: rmse 0.0002579027 max resid 0.0007659765
... Similar to previous best
Run 12 stress 0.1697659
Run 13 stress 0.1654879
Run 14 stress 0.1649254
Run 15 stress 0.1615097
Run 16 stress 0.1637516
Run 17 stress 0.1673051
Run 18 stress 0.1607993
Run 19 stress 0.1685889
Run 20 stress 0.1688124
*** Solution reachedscl <- 3 ## scaling == 3
colvec <- c("red2", "mediumblue")
plot(ord, type = "n")
with(speciesNPK, points(ord, display = "sites", col = colvec,
pch = 21, bg = colvec))
text(ord, display = "sites", cex = 0.8, col = "darkcyan")
with(speciesNPK, legend("topleft", levels(Microsite), bty = "n",
col = colvec, pch = 21, pt.bg = colvec))
adonis(speciesNPK[11:31]~speciesNPK$Microsite,method="jaccard")
Call:
adonis(formula = speciesNPK[11:31] ~ speciesNPK$Microsite, method = "jaccard")
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
speciesNPK$Microsite 1 0.8521 0.85212 3.9075 0.11874 0.001 ***
Residuals 29 6.3241 0.21807 0.88126
Total 30 7.1762 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(speciesNPK[11:31]~speciesNPK$NH4,method="jaccard")
Call:
adonis(formula = speciesNPK[11:31] ~ speciesNPK$NH4, method = "jaccard")
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
speciesNPK$NH4 1 0.3941 0.39415 1.6854 0.05492 0.068 .
Residuals 29 6.7821 0.23387 0.94508
Total 30 7.1762 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(speciesNPK[11:31]~speciesNPK$Clay_Depth,method="jaccard")
Call:
adonis(formula = speciesNPK[11:31] ~ speciesNPK$Clay_Depth, method = "jaccard")
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
speciesNPK$Clay_Depth 1 0.1646 0.16463 0.68089 0.02294 0.802
Residuals 29 7.0116 0.24178 0.97706
Total 30 7.1762 1.00000 Species do not show significant grouping by clay depth (p=.786) or ammonium (p=0.071) but do show significant grouping by microsite (p=0.001)
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] vegan_2.4-3 lattice_0.20-34 permute_0.9-4 MASS_7.3-47
[5] ggtern_2.2.0 ggplot2_2.2.1
loaded via a namespace (and not attached):
[1] latex2exp_0.4.0 Rcpp_0.12.10 DEoptimR_1.0-8
[4] plyr_1.8.4 compositions_1.40-1 tools_3.3.2
[7] boot_1.3-18 digest_0.6.12 nlme_3.1-129
[10] evaluate_0.10 tibble_1.3.3 gtable_0.2.0
[13] mgcv_1.8-16 rlang_0.1.1 Matrix_1.2-7.1
[16] yaml_2.1.14 parallel_3.3.2 proto_1.0.0
[19] gridExtra_2.2.1 cluster_2.0.5 stringr_1.2.0
[22] knitr_1.15.1 rprojroot_1.2 grid_3.3.2
[25] robustbase_0.92-7 bayesm_3.0-2 rmarkdown_1.4
[28] tensorA_0.36 magrittr_1.5 backports_1.0.5
[31] scales_0.4.1 htmltools_0.3.5 colorspace_1.3-2
[34] labeling_0.3 stringi_1.1.3 energy_1.7-0
[37] lazyeval_0.2.0 munsell_0.4.3 References
Brewer, J Stephen. 1999. “Effects of Fire, Competition and Soil Disturbances on Regeneration of a Carnivorous Plant (Drosera Capillaris).” The American Midland Naturalist 141 (1). BioOne: 28–42.
Craft, Christopher B, and Connie Chiang. 2002. “Forms and Amounts of Soil Nitrogen and Phosphorus Across a Longleaf Pine–depressional Wetland Landscape.” Soil Science Society of America Journal 66 (5). Soil Science Society: 1713–21.
Inglett, PW, KR Reddy, and R Corstanje. 2005. “Anaerobic Soils.” Encyclopedia of Soils in the Environment 1. Elsevier FL.
Olde Venterink, H, MJ Wassen, AWM Verkroost, and PC De Ruiter. 2003. “Species Richness–productivity Patterns Differ Between N-, P-, and K-Limited Wetlands.” Ecology 84 (8). Wiley Online Library: 2191–9.